Homework for Mr. Konstantin
Oksana Ayzsilnieks
March 29, 2017
library(DESeq2)
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
library(pheatmap)
library(outliers)
library(dplyr)
library(plotly)
library(grid)
library(futile.logger)
library(VennDiagram)
options(width=120)В ходе данной домашней работы понадобятся те же файлы, что и в лекции: “GSE89225_illumina_counts.csv”, “conditions.csv”, “human_mart.txt”. Для начала убедимся в том, что мы можем эти файлы прочитать. И посмотрим, что в них находится.
Неуместного аутлайера с условным обозначением treg_NBP_patient3 я вычислила, но немного позднее по ходу этого файла. Но, хотя Читатель узнает его имя позднее, я, раз знаю немного больше, удалю его из данных прямо сейчас сразу после считывания файла.
counts <- read.csv("GSE89225_Illumina_counts.csv", row.names=1)
conditions <- read.csv("conditions.csv", row.names=1)
mart <- read.table("human_mart.txt", sep="\t", header=1, check.names = F)
# а вот и его ненастоящее имя treg_NBP_patient3
counts <- select(counts, -treg_NBP_patient3)
conditions <- conditions[-12,]
# а это несколько строчек, чтобы взглянуть на струкутру датафреймов
print(counts[1:6, 1:2])## tconv_tumor_patient10 cd177negtreg_tumor_patient10
## ENSG00000000003 510 935
## ENSG00000000005 25 1
## ENSG00000000419 613 818
## ENSG00000000457 269 587
## ENSG00000000460 55 142
## ENSG00000000938 36 69dim(counts)## [1] 63677 33head(conditions)## tissue cells
## tconv_tumor_patient10 tissue: breast tumor cell type: Conventional CD4 T cells
## cd177negtreg_tumor_patient10 tissue: breast tumor cell type: Regulatory T cells
## cd177postreg_tumor_patient10 tissue: breast tumor cell type: Regulatory T cells
## tconv_tumor_patient1 tissue: breast tumor cell type: Conventional CD4 T cells
## treg_tumor_patient1 tissue: breast tumor cell type: Regulatory T cells
## tconv_NBP_patient2 tissue: NBP cell type: Conventional CD4 T cellsdim(conditions)## [1] 33 2head(mart)## Gene ID Associated Gene Name Gene type EntrezGene ID
## 1 ENSG00000252303 RNU6-280P snRNA NA
## 2 ENSG00000281771 Y_RNA misc_RNA NA
## 3 ENSG00000281256 RP11-222G7.2 lincRNA NA
## 4 ENSG00000283272 Clostridiales-1 sRNA NA
## 5 ENSG00000280864 RP11-654C22.2 lincRNA NA
## 6 ENSG00000280792 RP11-315F22.1 lincRNA NAdim(mart)## [1] 64101 4RNA-seq
Немного слов учителя об РНК-секвенировании:
Rna-seq steps:
- Изоляция РНК
- Rna selection / depletion
- Фрагментация
- Синтез кДНК
- Секвенирование
Rna selection / depletion:
- вся РНК
- тянем за поли(А)-хвосты (только мРНК)
- удаляем рибосомальную РНК (ribo-zero kit)
- таргетное секвенирование
Why Rna-seq?
- Не ограничены существующей сборкой и дизайном микрочипа
- Низкий фоновый сигнал
- Точность позволяет смотреть экспрессию отдельных изоформ
Sanity checks
###…и необходимой чистоте данных
Нужно всегда проверять длины библиотек и количество rRNA reads, которые оказались в библиотеке. Количество ридов можно проверять после выравнивания или после квантификации.
# Повторим сказанное на языке машины
proteinCoding <- mart[mart[, 3] == "protein_coding", ]
rRNA <- mart[mart[, 3] == "rRNA", ]
pcCounts <- counts[rownames(counts) %in% as.character(proteinCoding[, 1]), ]
rrnaCounts <- counts[rownames(counts) %in% as.character(rRNA[, 1]), ]
sampleCount <- ncol(counts)
toPlot <- data.frame(
sample=rep(colnames(counts), 3),
value=c(colSums(counts) - colSums(pcCounts) - colSums(rrnaCounts),
colSums(pcCounts),
colSums(rrnaCounts)),
type=c(rep("other", sampleCount),
rep("protein coding", sampleCount),
rep("rrna", sampleCount))
)
plot <- ggplot(data=toPlot, aes(x=sample, y=value, fill=type)) +
geom_bar(stat="identity") + theme_bw() +
theme(axis.text.x = element_text(angle=90, vjust=0.5))
plot
DESeq2
DESeq2 и правда хорош и, возможно, даже красив. Тут и дифференциальная экспрессия, и нормализации, и PCA-plots.
# исправилась, была неправа: вот теперь два DESeq-объекта
dds <- DESeqDataSetFromMatrix(countData = counts,
colData = conditions,
design = ~ tissue + cells)
dds <- dds[rowSums(counts(dds)) > 20, ]
dds <- DESeq(dds)
vst_dds <- vst(dds)
counts.norm <- assay(vst_dds)
# а вот и тот самый вышеупомянутый outlier, встречаем:
outlier(counts.norm)## tconv_tumor_patient10 cd177negtreg_tumor_patient10 cd177postreg_tumor_patient10 tconv_tumor_patient1
## 18.31254 19.06682 18.75782 19.50822
## treg_tumor_patient1 tconv_NBP_patient2 treg_NBP_patient2 tconv_tumor_patient2
## 19.42026 18.74975 18.97194 18.68064
## cd177negtreg_tumor_patient2 cd177postreg_tumor_patient2 tconv_NBP_patient3 tconv_tumor_patient3
## 18.86007 19.00433 18.76385 19.09414
## cd177negtreg_tumor_patient3 cd177postreg_tumor_patient3 tconv_NBP_patient4 treg_NBP_patient4
## 18.73277 19.78128 18.74858 19.24072
## tconv_tumor_patient4 cd177negtreg_tumor_patient4 cd177postreg_tumor_patient4 treg_tumor_patient5
## 18.51430 19.73984 19.74446 18.87674
## tconv_tumor_patient5 tconv_tumor_patient6 cd177negtreg_tumor_patient6 cd177postreg_tumor_patient6
## 18.84609 18.85645 18.68582 18.93730
## tconv_NBP_patient7 treg_NBP_patient7 tconv_NBP_patient8 treg_NBP_patient8
## 18.55074 18.03944 18.73227 18.73946
## treg_tumor_patient8 tconv_tumor_patient8 tconv_NBP_patient9 treg_NBP_patient9
## 18.43248 18.40272 18.90248 18.56628
## tconv_tumor_patient9
## 18.94650max(outlier(counts.norm)) # а это treg_NBP_patient3## [1] 19.78128# Вот такой еще
dds_cells <- DESeqDataSetFromMatrix(countData = counts,
colData = conditions,
design = ~ cells + tissue)
dds_cells <- dds_cells[rowSums(counts(dds_cells)) > 20, ]
dds_cells <- DESeq(dds_cells)
vst_dds_cells <- vst(dds_cells)
counts.norm_cells <- assay(vst_dds_cells)Поглядим на PCA без выброса?
pca_data_tissue <- prcomp(t(counts.norm))
percents_tissue <- pca_data_tissue$sdev^2 / sum(pca_data_tissue$sdev^2)
to_plot_tissue <- t(counts.norm) %*% pca_data_tissue$rotation
pca_data_cells <- prcomp(t(counts.norm_cells))
percents_cells <- pca_data_cells$sdev^2 / sum(pca_data_cells$sdev^2)
to_plot_cells <- t(counts.norm_cells) %*% pca_data_cells$rotation
gdata_tissue <- data.frame(
x=to_plot_tissue[, 1],
y=to_plot_tissue[, 2],
tissue=conditions[, 1],
name=rownames(conditions)
)
gdata_cells <- data.frame(
x=to_plot_cells[, 1],
y=to_plot_cells[, 2],
cells=conditions[, 2],
name=rownames(conditions)
)
# немного ненавязчивого интерактива от plotly
ggplotly(ggplot(data=gdata_tissue, aes(x=x, y=y, color=tissue, shape=tissue, text=name)) +
geom_point(size=3) + theme_bw() +
xlab(paste0("PC", 1, ": ", formatC(100 * percents_tissue[1], digits=4), "%")) +
ylab(paste0("PC", 2, ": ", formatC(100 * percents_tissue[2], digits=4), "%")))ggplotly(ggplot(data=gdata_cells, aes(x=x, y=y, color=cells, shape=cells, text=name)) +
geom_point(size=3) + theme_bw() +
xlab(paste0("PC", 1, ": ", formatC(100 * percents_cells[1], digits=4), "%")) +
ylab(paste0("PC", 2, ": ", formatC(100 * percents_cells[2], digits=4), "%")))ggplotly(plotPCA(vst_dds, intgroup=c("tissue", "cells"), ntop = 500) + theme_bw())Differential expression
“Давайте посмотрим, как выглядят результаты дифференциальной экспрессии и отсортируем их по статистике”.
res <- results(dds)
res <- res[order(res[,4]),]
res_cells <- results(dds_cells)
res_cells <- res_cells[order(res_cells[, 4]), ]Наконец исправила Volcano. А знаешь, что было не так? В объекте dds_cells на 135 строчке у меня была русская с!!!!!!!!!!!!!!!!!!!!
genes_tissue <- rownames(res)
genes_cells <- rownames(res_cells)
gdata_tissue <- data.frame(
log_fold_change=res$log2FoldChange,
p_adj=res$padj,
n = "Treg vs Tconv"
)
gdata_tissue <- na.omit(gdata_tissue)
gdata_cells <- data.frame(
log_fold_change=res_cells$log2FoldChange,
p_adj=res_cells$padj,
n = "Breast tumor vs Normal breast tissue"
)
gdata_cells <- na.omit(gdata_cells)
# я все это делаю еще и потому, что потом мне отдельно понадобятся датафреймы для клеток и тканей
# но для volcano, хотя я и планировала сначала использовать grid_arrange (почему и создавала два gdata для клеток и для тканей), в конечном итоге я использую слитый воедино датафрейм general_gdata. Потому как в последующем применю facet_grid.
general_gdata <- rbind(gdata_tissue, gdata_cells)
general_gdata <- na.omit(general_gdata)
# создадим порог значимости
general_gdata$threshold <- general_gdata$p_adj <= 0.01
general_gdata$significant <- sapply(general_gdata$threshold, function(x) ifelse(x == TRUE, "Significant", "Not Significant"))
general_gdata <- na.omit(general_gdata)
# Вы только посмотрите, какая красота!
ggplot(general_gdata,aes(x=log_fold_change, y=-log10(p_adj), col = significant)) +
facet_grid(. ~ n) +
geom_point(size=1) +
theme_bw() +
scale_colour_manual(values = c("black", "red")) +
labs(x = "Log fold change", y ="Adjusted p.value")+
geom_hline(yintercept=2, linetype = 2, size = 1, colour = 'red')
Поскольку я задержалась с домашкой, я решила, а может, еще что-нибудь почитать и как-нибудь поиграть с построением volcano? Прием, которым я воспользовалась ниже, направлен на более тонкое выявление дифференциальной экспрессии, основанное и на значении log2 fold change между двумя уровнями факторной переменной (уровнем экспрессии генов в разных клетках и в разных тканях) и на значении p-adj (p-val после поправки BH (Benjamini & Hochberg (1995)). Из полученного графика мы видим, что хотя различия в экспресии по показателю log2 fold change между некоторыми генами и выявляются, они остаются тем не менее незначимыми (черным цветом).
general_res <- rbind(res, res_cells)
p_adj <- c(general_res$padj)
log_fold_change <- c(general_res$log2FoldChange)
this_frame <- general_gdata[,-c(4,5)]
this_frame_green <- subset(this_frame, log_fold_change < -1 & p_adj < .01) # определим зеленую зону
this_frame_green <- cbind(this_frame_green, rep(1, nrow(this_frame_green)))
colnames(this_frame_green)[4] <- "Color"
this_frame_black <- subset(this_frame, (log_fold_change >= -1 & log_fold_change <= 1) | p_adj >= .01) # определим плохую черную зону
this_frame_black <- cbind(this_frame_black, rep(2, nrow(this_frame_black)))
colnames(this_frame_black)[4] <- "Color"
this_frame_red <- subset(this_frame, log_fold_change > 1 & p_adj < .01) # определим красную зону
this_frame_red <- cbind(this_frame_red, rep(3, nrow(this_frame_red)))
colnames(this_frame_red)[4] <- "Color"
frame_total <- rbind(this_frame_green, this_frame_black, this_frame_red)
frame_total$Color <- as.factor(frame_total$Color)
genes <- rownames(general_res)
##Сконструируем объект volcano
my_another_volcano_hello <- ggplot(data = frame_total, aes(x = log_fold_change, y = -log10(p_adj),col= Color)) +
geom_point(alpha = 0.5, size = 1) + theme_bw() + theme(legend.position = "none") +
geom_hline(yintercept = 2, colour = "black", linetype = 2, size = 1) +
scale_color_manual(values = c("green", "black", "red")) +
facet_grid(. ~ n) +
labs(y="Adjusted p.value", x="Log fold change")
ggplotly(my_another_volcano_hello)А теперь взглянем на горячую карту. Горячая карта покажет горячие точки. Вдруг эти точки о чем-то говорят. (Шучу, - конечно же, они говорят)
# процедуры извлечения генов из пасвея
kkeys <- keys(org.Hs.eg.db, keytype="ENSEMBL")
goAnno <- AnnotationDbi::select(org.Hs.eg.db, keys=kkeys,
keytype="ENSEMBL", columns=c("GOALL", "ONTOLOGYALL", "SYMBOL"))
goAnno <- tbl_df(goAnno)
goAnno <- filter(goAnno, GOALL=="GO:0007159")
# or you can pick ENTREZ, or SYMBOL, or whatever you want
genesToVisualise <- goAnno$ENSEMBL
# несколько шагов перед тем, как предложить что-то изящной функции
to_visualise_zzzz <- counts.norm[rownames(res), order(conditions[, 2])]
subsetting <- rownames(subset(res, rownames(res) %in% genesToVisualise))
nearly_want_to_see <- to_visualise_zzzz[subsetting,]
exactly_want_to_see <- t(apply(nearly_want_to_see, 1, function(r) {(r - min(r)) / (max(r) - min(r))}))
# О, как она изящна
pheatmap(exactly_want_to_see,
show_rownames = F, cluster_rows = F,
cluster_cols=F,
annotation_col = conditions,
main = "GO:0007159: leukocyte cell-cell adhesion")
Построим же наконец диаграмму Венна
for_set_1 <- rownames(res[which(res$padj < 0.01),])
for_set_2 <- rownames(res_cells[which(res_cells$padj < 0.01),])
wanted_sected <- intersect(for_set_1, for_set_2)
length(wanted_sected) # 84## [1] 84venn_plot <- draw.pairwise.venn(length(for_set_1), length(for_set_2), length(wanted_sected), category = c("Tregs vs Tconv", "Tissues"), fill = c("blue", "pink"), ext.line.lty = "dashed")